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Abstract 


This  report  presents  a  fast  hybrid  volume-surface  integral  equation  approach  for  the  computation 
of  electromagnetic  scattering  from  composite  left-handed  materials  (LHM)  such  as  split-ring  res¬ 
onators  (SRR)  with  wires.  The  volume  electric  field  integral  equation  (EFIE)  is  applied  to  the 
dielectric  region  of  this  LHM,  and  the  surface  electric  field  integral  equation  is  applied  on  the 
conducting  surface.  The  method  of  moments  (MoM)  is  used  to  discretize  the  integral  equation 
into  a  matrix  solution  and  adaptive  integral  method  (AIM)  is  employed  to  reduce  the  memory 
requirement  and  CPU  time  for  the  matrix  solution.  The  present  approach  is  sufficiently  versatile 
in  handling  scattering  problems  of  such  composite  LHMs,  due  to  the  combination  of  volume  and 
surface  electric  field  integral  equations.  This  algorithm  is  first  verified  in  a  scattering  problem  of 
an  ordinary  composite  right-handed  material  structure.  Then  we  use  the  algorithm  to  find  the 
resonance  frequency  of  a  few  SRRs.  In  order  to  demonstrate  the  accuracy  and  efficiency  of  this 
technique,  we  present  the  numerical  results  of  calculating  the  radar  cross  section  (RCS)  of  such  an 
SRR  slab.  Then  based  on  previous  knowledge  of  such  structure,  we  designed  a  new  left-handed 
structure  on  microstrip  line  with  capacitor  gaps.  Simulation  results  show  the  backward- wave  char¬ 
acteristic  of  this  structure.  And  then  we  apply  this  structure  in  antenna  design  to  improve  the 
performance  of  an  ordinary  microstrip-feed  patch  antenna.  The  antenna’s  bandwidth  is  greatly 
broadened  while  the  VSWR  remains  very  low  and  radiation  efficiency  is  high.  In  the  end,  some 
future  research  proposals  focus  on  interesting  and  practical  topics  are  also  included  in  the  report. 
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Chapter  1 


Introduction 


Recently,  left-handed  material(LHM),  which  is  also  known  as  negative  index  media  (NIM),  has 
attracted  a  great  deal  of  attention.  These  materials  exhibit  a  negative  refraction  index  n,  and 
therefore  can  have  very  interesting  properties,  e.g.,  the  phase  and  group  velocities  of  electromag¬ 
netic  waves  in  this  material  are  opposite.  Since  the  E  and  H  fields  and  the  wave  vector  k  of  a 
propagating  plane  EM  wave  form  a  left-handed  system  in  these  materials,  Veselago  [1]  introduced 
the  term  “left-handed  materials”  in  his  seminal  work  published  in  1968.  The  interest  in  Veselago’s 
work  was  renewed  since  Pendry  et  al.  [2]  proposed  an  artificial  material  consisting  of  the  so  called 
split-ring  resonators  (SRRs)  as  shown  in  Figure  1.2,  which  exhibit  a  band  of  negative  //  values  in 
spite  of  being  made  of  non-magnetic  materials  and  wires  which  provide  the  negative  e  behavior. 
Based  on  Pendry’s  suggestions  and  targeting  the  original  idea  of  Veselago,  Smith  et  al.  demon¬ 
strated  in  2000  the  realization  of  the  first  left-handed  material  which  consisted  of  an  array  of  split 
ring  resonators  (SRRs)  and  thin  wires  in  alternating  layers  [3].  To  date,  the  range  of  imaginable 
applications  of  left-handed  materials  extends  to  the  field  of  antenna  design,  vehicle  coatings  for 
altering  radar  cross  section  properties,  and  lenses  [4]. 

Since  the  original  microwave  experiment  of  Smith  et  al.,  several  LHMs  were  fabricated  [5]  [6]  that 
exhibited  a  pass  band  in  which  it  was  assumed  that  e  and  /j  are  both  negative.  This  assumption 
was  based  on  transmission  measurements  of  the  wires  alone,  the  SRR  alone,  and  the  LHMs.  The 
occurrence  of  a  LHM  transmission  peak  within  the  stop  bands  of  the  SRR  and  wire  structures 
was  taken  as  evidence  for  the  appearance  of  LHM  behavior.  Further  support  to  this  interpretation 
was  provided  by  the  demonstration  that  such  LHMs  exhibit  negative  refraction  of  EM  waves  [7]  [8]. 
Moreover,  there  is  a  significant  amount  of  numerical  work  [9]  [10]  [11]  [12]  in  which  the  transmission 
and  reflection  data  are  calculated  for  a  finite  length  of  metamaterial.  A  retrieval  procedure  can 
then  be  applied  to  obtain  the  metanraterial  parameters  £  and  //  under  the  assumption  that  it  can 
be  treated  as  homogeneous.  This  procedure  was  applied  in  [13]  and  confirmed  that  a  medium 
composed  of  SRRs  and  wires  could  indeed  be  characterized  by  effective  £  and  /i  whose  real  parts 
were  both  negative  over  a  finite  frequency  band,  as  was  the  real  part  of  the  refractive  index  n. 
However,  it  was  recently  shown  [14]  that  the  SRRs  exhibit  a  resonant  electric  response  in  addition 
to  their  resonant  magnetic  response,  and  therefore,  make  the  understanding  and  the  fabrication  of 
true  LHMs  more  challenging. 
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Poynting’s  vector  S 
Vector  k  1 


Source  of  radiation 


Figure  1.1:  Poyntings  vector  and  the  wave  vector  for  wave  propagations  in  right-handed  and  left- 
handed  materials. 


Figure  1.2:  Round  and  square  split-ring  resonators 
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To  design  such  artificial  materials  for  different  specific  engineering  applications,  simulation  tools 
are  essential  for  saving  efforts  and  costs.  The  numerical  methods  for  characterizing  these  struc¬ 
tures  are  still  those  fundamental  numerical  simulation  techniques,  i.e. ,  the  differential-equation 
(DE)  method  [15]  [16],  the  integral-equation  (IE)  method  [17] [18],  and  the  equivalent-circuit  models 
[19]  [20].  Among  the  three  methods  usually  utilized  for  designing  artificial  synthetic  materials,  the 
electromagnetic  solver  developed  based  on  the  method  of  moments  (MoM)  for  solving  the  hybrid 
volume-surface  integral  equation  (VSIE)[21]  is  rigorous  and  widely  adopted  in  various  electromag¬ 
netic  simulations.  It  is  certainly  a  very  good  option  for  analyzing  electromagnetic  characteristics 
of  metanraterials.  The  method  of  moment  is  one  of  the  most  popular  numerical  method  for  solving 
boundary-value  problems  in  electromagnetic. 

However,  the  MoM  will  lead  to  a  full  matrix  whose  numerical  iterative  solution  requires  0{NnerN 2) 
operations  and  0(N2)  memory  to  store  the  matrix  elements,  where  N  and  Niter  denote  the  num¬ 
ber  of  unknows  and  number  of  iterations,  respectively  [22].  The  computational  cost  and  memory 
requirement  will  increase  much  faster  with  the  larger  number  of  unknowns.  On  one  hand,  the 
operating  frequency  of  electromagnetic  waves  in  a  medium  of  metanraterials  falls  usually  within 
the  microwave  frequency  band  or  even  higher.  On  the  other  hand,  the  spacings  and  dimensions 
of  conducting  inclusions  that  compose  the  artificial  media  are  required  to  be  much  smaller  than 
the  wave-length  in  order  to  synthesize  a  practical  continuous  material  so  that  the  number  of  in¬ 
clusions  is  numerous  and  the  resultant  total  number  of  unknowns  is  very  large,  even  though  the 
dimensions  of  such  synthesized  metanraterials  are  electrically  small  in  this  frequency  band.  To 
achieve  this  objective,  a  fast  algorithm  has  to  be  integrated  into  the  MoM  procedure  to  accelerate 
the  matrix- vector  multiplications. 

In  Chapter  2,  we  discussed  some  special  characteristic  of  split-ring  resonators  as  well  as  the 
adaptive  integral  method  (AIM).  Here  we  AIM  into  the  scattering  problem  of  LHM.  It  was  a  fast 
integral-equation  solver  applicable  to  large-scale  electromagnetic  scattering  problems  which  was 
described  by  Bleszynski  [23]  [24]  for  the  first  time  in  1996.  In  this  approach,  the  field  generated  by 
a  given  current  distribution  is  decomposed  into  near  and  far  field  components.  The  near  field  is 
computed  using  the  conventional  method  of  moments  technique  with  the  Galerkin  discretization. 
The  far  field  is  calculated  by  approximating  the  original  current  distribution  by  an  equivalent 
current  distribution  on  a  regular  Cartesian  grid,  such  that  the  two  currents  have  identical  multipole 
moments  up  to  a  required  order  m.  As  the  result  of  this  discretization,  the  original  full  impedance 
matrix  is  decomposed  into  a  sum  of  a  sparse  matrix  (corresponding  to  the  near  field  component) 
and  a  product  of  sparse  and  Toeplits  matrices  [25]  (corresponding  to  the  far  field  component). 
Because  of  the  convolution  nature  of  the  Toeplits  kernel,  the  field  generated  by  the  equivalent 
current  distribution  can  be  then  obtained  by  means  of  discrete  Fast  Fourier  transforms. 

In  Chapter  3,  we  show  the  numerical  result  of  our  simulations  based  on  this  AIM  algorithm. 
First,  we  design  the  size  and  shape  of  our  structure.  Then  with  AIM,  we  calculated  the  radar 
cross  section  (RCS)  of  a  few  SRR  inclusions  at  an  estimated  frequency  range  to  find  the  resonance 
frequency  of  this  particular  structure.  At  the  vicinity  of  this  frequency,  we  examine  the  RCS 
of  a  SRR  slab  with  different  kinds  of  plane-wave  input.  From  the  result,  we  find  that  when 
the  input  plane  wave  satisfy  the  magnetic  resonance  condition  for  SRR  [14],  the  power  can  be 
transferred  which  satisfy  the  property  of  double-negative.  However,  when  the  input  doesn’t  satisfy 
the  resonance  condition,  the  scattered  electric  field  is  very  weak  since  the  whole  slab  becomes  no 
more  than  a  ordinary  scatterer.  Then  on  one  hand,  this  phenomenon  obtained  re-proved  the  LH 
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property  of  this  SRR  slab,  on  the  other  hand,  it  confirms  that  our  AIM  algorithm  is  an  efficient 
and  accurate  solver  to  analyze  composite  LHM. 

In  Chapter  4,  we  present  some  newly  designed  left-handed  structures  and  make  necessary  im¬ 
provement  to  broaden  the  bandwidth.  The  new  structure  is  applied  in  the  fabrication  of  patch 
antenna  to  enhance  the  performance.  Good  agreements  are  achieved  between  simulation  results 
and  practical  measurements.  The  antenna  performs  well  with  high  gain,  low  loss  and  wide  band. 

In  Chapter  5,  future  research  topics  and  possible  method  is  proposed. 


Chapter  2 


Split-Ring  Resonators  and  Adaptive 
Integral  Method 


2.1  Split-Ring  Resonators 


Currently,  the  metamaterials  under  study  are  composed  of  a  periodic  arrangement  of  metallic 
printed  lines  (or  rods)  that  exhibit  a  negative  permittivity,  and  split-ring  resonators  (SRRs)  that 
exhibit  a  negative  permeability.  To  understand  how  these  nretanraterials  function,  it  is  best  to 
follow  the  historical  steps  that  lead  to  the  current  geometry.  First,  it  was  shown  that  a  3-D  array 
of  thin,  continuous  (very  long)  metal  wires  exhibit  a  resonant  frequency  response  similar  to  that  of 
a  plasma  medium  [26]  [27] .  At  frequencies  below  the  resonant  frequency,  the  permittivity  becomes 
negative,  the  wave  vector  is  imaginary,  and  there  is  no  transmission.  Above  the  resonant  frequency, 
the  permittivity  is  real  and  transmission  occurs.  Second,  it  was  demonstrated  that  negative  per¬ 
meability  could  be  achieved  using  an  array  of  SRRs  [2].  The  SRRs  act  like  magnetic  dipoles,  with 
a  resonant  response  resulting  from  internal  inductance  and  capacitance.  Close  to  resonance,  the 
SRRs  produce  strong  magnetic  fields  whose  directions  oppose  the  incident  magnetic  field,  resulting 
in  a  negative  effective  permeability.  Finally,  once  the  negative  permeability  structures  had  been 
discovered,  finding  ways  of  combining  them  with  negative  permittivity  structures  was  the  next 
step  in  research,  leading  to  the  well  known  geometry  shown  in  [28].  However,  the  demonstration 
of  the  LH  properties  of  a  metanraterial  is  not  as  straightforward  as  it  seems,  and  showing  that 
permittivity  and  permeability  are  separately  negative  over  a  similar  frequency  band  does  not  im¬ 
ply  that  they  remain  so  when  the  rods  and  SRRs  are  combined  into  a  unique  geometry.  One  of 
the  first  observations  of  the  LH  nature  of  these  metanraterial  structures  was  inferred  by  observing 
transmission  and  reflection  coefficients  over  a  wide  frequency  band.  For  the  separate  geometries 
(rods  only  or  SRRs  only),  as  shown  in  Figure  2.3  no  transmission  occurs  because  of  an  imagi¬ 
nary  wave  number.  In  the  frequency  band  where  the  two  parameters  are  simultaneously  negative, 
transmission  occurs  when  the  combined  geometry  is  considered.  This  band-pass  phenomenon  was 
indeed  observed  in  these  nretanraterials  [28].  However,  as  mentioned  before,  this  transmission  is 
not  sufficient  to  conclude  that  the  material  is  exhibiting  LH  behavior,  since  the  interaction  between 
the  SRRs  and  rods  could  result  in  a  positive  permittivity  and  permeability,  which  would  also  yield 
transmission.  To  solve  this  problem  Moss  et  al.  did  a  more  rigorous  test  to  demonstrate  the  LH 
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Figure  2.1:  Split-rings  resonators  with  rods 


Figure  2.2:  Split-ring  resonators  with  wires 
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Plasma  and  line 
wire  structures 


e  <  0,  p  >  0 
No  transmission 


Microstructured 
magnets  and 
split  rings 


e  >  0,  p  <  0 

No  transmission 


Figure  2.3:  Different  combination  of  e  and  p 

nature  of  a  material  [29]  which  show  that  the  combination  of  SRRs  and  rods  does  yield  a  pass 
band  of  double  negative.  So,  base  on  all  these  research  result,  we  use  SRR  structure  in  further 
modeling  and  simulation. 


2.2  Method  of  Moments 


To  design  and  compute  the  EM  field  of  such  artificial  materials,  simulation  tools  are  essential 
for  saving  efforts  and  costs.  The  numerical  methods  for  characterizing  these  structures  are  still 
those  fundamental  numerical  simulation  techniques,  i.e. ,  the  differential-equation  (DE)  method 
[15] [16], the  integral-equation  (IE)  method  [17]  [18],  and  the  equivalent-circuit  models  [19]  [20].  Among 
the  three  methods  usually  utilized  for  designing  artificial  synthetic  materials,  the  electromagnetic 
solver  developed  based  on  the  method  of  moments  (MoM)  for  solving  the  hybrid  volume-surface 
integral  equation  (VSIE)[21]  is  rigorous  and  widely  adopted  in  various  electromagnetic  simulations. 

The  method  of  moments  (MoM)  is  a  general  procedure  for  solving  linear  equations.  Consider 
the  inhomogeneous  equation 


L(f)  =  9  (2.1) 

where  L  is  a  linear  operator,  g  is  known,  and  /  is  to  be  determined.  Let  /  be  expanded  in  a  series 
of  functions  /i,  /2,  /3,  ■  ■  •  in  the  domain  of  L,  as 

f  =  Yla-fn  (2-2) 

n 

where  the  an  are  constants.  We  shall  call  the  fn  expansion  functions  or  basis  functions.  For  exact 
solution,  (2.2)  is  usually  an  infinite  summation  and  the  fn  form  a  complete  set  of  basis  functions. 
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For  approximate  solution,  (2.2)  is  a  finite  summation.  Substituting  (2.2)  in  (2.1),  and  using  the 
linearity  of  L,  we  have 

]T  anL(fn)  =  g  (2.3) 

n 

It  is  assumed  that  a  suitable  inner  product  <  f,g  >  has  been  determined  for  the  problem.  Now 
define  a  set  of  weighting  functions,  or  testing  functions,  wi,w2,ws,  ...  in  the  range  of  L  and  take 
the  inner  product  of  (2.3)  with  each  wm.  The  result  is 

^  ^  Q^n  ^  tr /n ■  /.  f a  > — <  wm,g  >  (2-4) 

n 

m  =  1,2,3,...  This  set  of  equations  can  be  written  in  matrix  form  as 


where 


[^mn]  [c^n]  —  [dm] 


<wi,Lfi>  <w1,Lf2> 
<w2,Lfi>  <  w2l  Lf2  > 


aq 

<  W\. 

9>1 

.C^n]  — 

OL2 

\lhn\  — 

<  w2, 

1 

A 

Ol 

(2.5) 


(2.6) 


(2.7) 


If  the  matrix  [/]  is  nonsingular  it  inverse  [l  x]  exists.  The  an  are  then  given  by 

[an\  =  [l nm][Srn. ) 


(2.8) 


and  the  solution  for  /  is  given  by  (2.2).  For  concise  expression  of  this  result,  define  the  matrix  of 
functions 


[fn]  =  [fl  h  h  ■  ■]  (2.9) 

and  write 

/=[/n][U]M  (2-10) 

This  solution  may  be  exact  or  approximate,  depending  upon  the  choice  of  the  fn  and  wn.  The 
particular  choice  wn  =  fn  is  known  as  the  Galerkin’s  method. 

If  the  matrix  [/]  is  of  infinite  order,  it  can  be  inverted  only  in  special  cases,  for  example,  if  it  is 
diagonal.  The  classical  eigenfunction  method  leads  to  a  diagonal  matrix,  and  can  be  thought  of  as 
a  special  case  of  the  method  of  moments.  If  the  sets  fn  and  wn  are  finite,  the  matrix  is  of  finite 
order,  and  can  be  inverted  by  known  method. 

One  of  the  main  tasks  in  may  particular  problem  is  the  choice  of  the  fn  and  wn.  The  fn  should 
be  linearly  independent  and  chosen  so  that  some  superposition  can  be  approximated  /  reasonably 
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well.  The  wn  should  also  be  linearly  chosen  so  that  the  products  <  wn,g  >  depend  on  relatively 
independent  properties  of  g.  Some  additional  factors  which  affect  the  choice  of  fn  and  wn  are  the 
accuracy  of  solution  desired,  the  ease  of  evaluation  of  the  matrix  elements,  the  size  of  the  matrix 
that  can  be  inverted,  and  the  realization  of  a  well-conditioned  matrix  [/] 

However,  the  MoM  will  lead  to  a  full  matrix.  The  matrix  equation  can  be  solved  using  a 
conventional  direct  solver  or  an  iterative  solver.  A  direct  solver  requires  0(N 3)  operations  to  solve 
the  equation  while  an  iterative  solver  requires  0{NnerN 2)  operations.  Both  solvers  require  0(N2) 
memory  storage  to  store  the  matrix  elements  where  N  and  Nner  denote  the  number  of  unknows 
and  number  of  iterations,  respectively  [22].  These  stringent  requirements  have  prevented  the  use 
of  the  conventional  MoM  to  solve  problems  with  a  large  number  of  unknowns.  The  computational 
cost  and  memory  requirement  will  increase  much  faster  with  the  larger  number  of  unknowns.  On 
one  hand,  the  operating  frequency  of  electromagnetic  waves  in  a  medium  of  metamaterials  falls 
usually  within  the  microwave  frequency  band  or  even  higher.  On  the  other  hand,  the  spacings  and 
dimensions  of  conducting  inclusions  that  compose  the  artificial  media  are  required  to  be  smaller 
than  the  wave-length  in  order  to  synthesize  a  practical  continuous  material  so  that  the  number 
of  inclusions  is  numerous  and  the  resultant  total  number  of  unknowns  is  very  large,  even  though 
the  dimensions  of  such  synthesized  metamaterials  are  electrically  small  in  this  frequency  band.  To 
achieve  this  objective,  a  fast  algorithm  has  to  be  integrated  into  the  MoM  procedure  to  accelerate 
the  matrix- vector  multiplications. 


2.3  Adaptive  Integral  Method 

2.3.1  the  Advantages  of  Adaptive  Integral  Method  (AIM) 

There  are  three  major  fast  algorithms  recently  developed  to  overcome  the  problem  of  large  memory 
requirement  and  computational  complexity  caused  by  MoM. 

1)  Fast  multipole  method  (FMM)  [30,  31]  and  its  extension,  multilevel  fast  multipole  algorithm 
(MLFMA)  [32],  exhibit  computational  complexity  on  the  order  of  0(N 1"5)  and  0(N  log  N),  re¬ 
spectively.  However,  a  large  value  of  N  is  required  to  gain  advantage  over  the  MoM,  and  certainly 
this  is  satisfied  mostly  when  the  FMM  or  MLFMA  is  employed. 

2)  Conjugate  gradient  fast  Fourier  transform  (CG-FFT)  as  a  grid-based  fast  algorithm  has 
been  also  used  widely  to  solve  the  scattering  problem  formulated  by  surface  integral  equation 
(SIE)  and  volume  integral  equation  (VIE)  [33].  The  attracting  feature  of  CG-FFT  is  the  memory 
requirement  and  computational  complexity  are  of  O(NlogN)  and  O(N),  respectively.  However, 
the  major  drawback  of  this  method  is  the  staircase  approximation  due  to  the  usage  of  uniform 
rectangular  grid  to  model  arbitrarily  shaped  object,  and  therefore  the  resulted  staircase  error. 

3)  The  adaptive  integral  method  (AIM)  has  been  proposed  to  overcome  the  shortcomings  of 
CG-FFT.  It  was  initially  proposed  to  solve  scattering  problem  of  conducting  objects  [23,  24],  The 
AIM  and  its  counterpart,  precorrected  FFT  (P-FFT),  were  also  used  to  solve  the  other  problem 
formulated  using  the  surface  integral  equation  (SIE)  [34] .  Recently,  Zhang  and  Liu  have  used  the 
AIM  to  solve  the  scattering  problem  formulated  using  VIE  [35].  The  computational  complexity 
for  AIM  is  0(N1'5  log  N)  and  0(N  log  N)  for  surface  and  volume  scatterers,  respectively. 
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Later  in  this  report,  this  AIM  will  be  applied  for  computing  the  electromagnetic  scattering  by 
the  SRR  structures  to  examine  their  LH  property. 


2.3.2  Introduction  of  Hybrid  Volume-Surface  Integral  Equation 

The  scattering  problem  of  mixed  dielectric  and  conducting  objects  can  be  formulated  using  SIE 
or  the  combination  VIE  and  SIE.  The  scattering  by  conducting  objects  is  usually  characterized 
by  using  the  SIE,  but  the  choice  of  using  SIE  or  VIE  for  dielectric  object  depends  on  its  inhomo¬ 
geneity.  For  an  object  with  arbitrary  inhonrogeneity,  solving  the  scattering  problem  using  VIE  will 
be  more  advantageous  than  SIE,  as  the  latter  requires  the  integral  equations  to  be  formulated  in 
every  dielectric  regions.  However,  if  the  object  consists  of  only  piecewise  homogeneous  dielectric 
materials,  then  SIE  formulation  is  generally  preferred.  The  formulation  of  mixed  dielectric  and 
conducting  object  problem  using  only  SIE  can  be  found  in  the  literature  [36]  [37].  Another  for¬ 
mulation  that  involves  both  the  VIE  and  SIE  was  made  to  describe  the  scattering  by  dielectric 
and  conducting  bodies  and  then  to  combine  them  together  to  form  hybrid  volume-surface  integral 
equation  (VSIE)  [38] .  We  found  this  method  is  useful  to  SRR  structure  since  VIE  and  SIE  can  be 
applied  to  describe  the  EM  field  in  the  dielectric  region  and  conducting  surface  separately. 

And,  for  closed  conducting  body,  the  combined  field  integral  equation  (CFIE)  can  be  used  as  it 
eliminates  interior  resonance  problem  and  it  normally  yields  a  better  conditioned  matrix  equation 
as  compared  to  electric  field  integral  equation  (EFIE)  or  magnetic  field  integral  equation  (MFIE). 
However,  for  the  open  conducting  body  like  SRR,  the  MFIE  can  not  be  used.  So  here  we  use  EFIE 
only  for  SIE. 


2.3.3  Formulations  For  Hybrid  Volume-Surface  Integral  Equation 

In  this  section,  we  discuss  the  formulations  for  VSIE. 

First,  we  consider  a  general  case  of  an  arbitrarily  shaped  three-dimensional  (3-D)  scatterer, 
which  consists  of  homogeneous  or  inhomogeneous  dielectric  material  and  conducting  body  as  shown 
in  Figure  2.4.  The  object  is  embedded  in  an  isotropic  homogeneous  background  medium  with 
permittivity  £b  and  permeability  pb-  The  scatterer  is  illuminated  by  an  incident  wave  Emc,  which 
is  excited  by  impressed  sources  in  the  background  media.  The  dielectric  region  V  is  assumed  to 
have  the  permeability  p  =  pb  and  complex  dielectric  constant  e  =  e(r)  —  ja(r)/uj  where  e(r)  and 
cr  (?')  are  permittivity  and  conductivity,  respectively,  at  r.  Inside  the  dielectric  region  V  and  on  the 
surface  of  conducting  body  S,  the  incident  wave  induces  volume  current  J y  and  surface  current 
Js ■  The  scattered  electric  field  Esca  and  magnetic  field  Hsca  produced  by  the  induced  volume 
and  surface  currents  are  given  by 

E%a(r)  =  -jkbribAn(r)-V$n(r),  Q  =  S  or  V  (2.11) 

H'Qa(r)  =  V  x  An(r),  Q  =  S  or  V  (2.12) 

where  the  magnetic  vector  potential  is  defined  as: 


Afi(r)  =  fnJn(r')g(r,r')dr',  H  =  S  or  V 


(2.13) 
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Figure  2.4:  Geometry  of  scatterer  consisting  of  dielectric  material  and  conducting  body  embedded 
in  an  isotropic  homogeneous  medium 

and  the  electric  scalar  potential  is  defined  as: 

$n(r)  =  -jfcfnV'-Jn(r')g(r,r'dr'),  Q  =  S  or  V  (2.14) 

—jkfj \V  —  V'\  / 

with  g(r,r')  =  e47r|r_r/|  >  Vb  =  wf-A  and  kb  as  the  wavenumber  of  background  medium 

In  the  dielectric  region  V,  by  taking  the  scattered  field  from  the  surface  current  into  considera¬ 
tion,  the  total  electric  field  becomes 

E(r)  =  Einc(r)  +  EsJa(r)  +  Essca{r)  (2.15) 

and  the  induced  volume  current  is  defined  as 

Jv  =  ju(i  -  eb)E  (2.16) 


Similarly,  the  scattered  field  by  induced  volume  current  also  contributes  to  the  total  field  on  S. 
Since  the  tangential  components  of  total  electric  field  vanish  on  conducting  surface,  we  get 

h  x  Einc(r)  =  -h  x  EsJa(r)  -  h  x  Essca(r)  (2.17) 

Equations  2.15  and  2.17  are  known  as  the  EFIE  as  the  formulations  involve  only  electric  field.  For  a 
closed  conducting  surface,  the  MFIE  can  be  obtained  by  considering  the  tangential  components  of 
the  total  magnetic  field  on  the  conducting  surface,  which  are  equal  to  the  induced  surface  current 
components.  Since  we’ll  not  use  MFIE  here,  we  do  not  list  the  equations  here. 


2.3.4  RWG  and  SWG  basis  function 

As  we  mentioned  in  the  section  of  MoM,  we  must  have  basis  functions  to  expand  the  current 
on  conducting  surface  and  dielectric  volume.  These  basis  functions  must  be  associated  with  the 
surface  and  volume  element  we  get  after  discretization.  Here  in  our  problem,  the  volume  of  dielectric 
material  and  surface  of  conducting  body  are  discretized  into  tetrahedral  elements  and  triangular 
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patches,  respectively.  These  elements  are  used  because  of  their  flexibility  to  model  arbitrarily 
shaped  3-D  object.  Then  the  surface  current  is  expanded  using  vector  basis  function  fs ,  which  has 
constant-normal  and  linear-tangential  (CN/LT)  components  along  the  edge 

NS 

JS  =  J2Jnfn  (2-18) 

n=l 

As  for  the  dielectric  region,  the  electric  flux  D  is  expanded  using  appropriate  volume  CN /LT  basis 
functions  fy  as  follows: 

Nv 

D  =  i(r)E  =  J2lnfn  (2-19) 

n=  1 


by  use  of  equation  2.16,  the  induced  volume  current  is  expressed  as 


Nv  \ 

T  •  ~  £°  TV 

j V  =  JU  V V-— 

1  e(r) 

n=  1 


Nv 


Jnfn  =  Y  K(r)In  f 


n=  1 


(2.20) 


where  n(r)  = 


_  e(r)-e o 
£-(r) 


is  the  contrast  ratio  and  i(r)  is  the  permeability  of  a  tetrahedron  element. 


In  this  report,  the  CN/LT  basis  functions  that  have  been  used  are  SWG  basis  function  [39]  and 
RWG  basis  function  [40]  for  the  induced  volume  current  and  surface  current,  respectively.  The  SWG 
basis  function  is  defined  on  a  pair  of  tetrahedral  elements  while  the  RWG  basis  function  is  defined 
on  a  pair  of  triangular  patches.  The  properties  of  SWG  and  RWG  basis  functions  for  example, 
the  continuity  of  the  electric  flux  density  normal  to  the  interior  face  (SWG)  and  the  continuity  of 
the  surface  current  normal  to  the  interior  edge  (RWG)  make  them  suitable  to  efficiently  discretize 
the  previous  integral  equations.  The  formulas  for  SWG  and  RWG  basis  function  as  well  as  the 
testing  procedure  we  used  in  our  simulation  program  will  not  be  list  here  since  the  these  complicated 
formulas  are  not  our  main  point.  However,  they  can  be  found  in  the  Appendix.  It  is  also  noted  that 
during  discretization  of  the  object,  the  triangular  patches  are  generated  such  that  it  is  coincided 
with  the  surface  of  tetrahedron.  Using  the  volume  basis  function  to  test  (2.15)  and  the  surface 
basis  function  to  test  (2.17),  the  hybrid  integral  equations  will  be  converted  into  a  matrix  equation 
system  as: 


zss 


[  Iv  1 

'  Ev  ' 

Is 

Es 

(2.21) 


where  the  vectors  /'  and  Is  represent  the  coefficients  of  volume  current  and  surface  current 
respectively.  The  excitation  vector  can  be  computed  using 


yv  — 


f  /m  ’  Einc(r')dr' 

JVm 


Ebm  =  fbm-  Emc(r')dr' 

JSm 


(2.22) 

(2.23) 
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while  the  elements  of  the  block  matrices  in  (2.21)  can  be  computed  using  the  formulas  shown 
below: 


7vv  = 

^ mn 


7VS  = 

^ mn 


7SV  = 

^ mn 


7SS  — 

^ mn 


fV  nV 
—2—2- dr 


JVm  E 

- wkbr]b  i 


>Vm 


V 


m  i  /  '■■/!.  !l(  f-  r')dr'  +  T9/  V'  •  (k/£)  g[r,  r’) 

Kb  Jvn 


L  JVn 


jkbhb  /  fm 

■Jvm 

-ukbi]b  /  fSm 
■JSm 

jhdb  /  fm 

JSm 


V 


USn 


fn9{r,  r')dr  +  p  V  ■  f„g{r,  r')dr 


dr 


dr  (2.24a) 
(2.24b) 


/  Kfn9(r,  r')dr'  +  5  /  V'  •  (k/£ )  g{r,  r')dr' 

UVn  kb  JVn 


V 


/  fn9(r>  r')dr'  +  V2  v'  •  fn9{r,  r')dr' 
>Sn  H  Jsn 


dr 


dr  (2.24c) 
(2.24d) 


k  is  a  constant  value  within  the  tetrahedron  and  can  be  taken  out  from  the  integration. 
Besides,  we  can  also  write 

V-(k/J0=kV-/£  +  V«-/£ 

and  the  second  term  given  by  2.25  as: 

fv  _  \  dn  ( Kn  ~  k~),  on  the  common  facet 
K  n  1  0,  elsewhere 


(2.25) 


(2.26) 


2.3.5  Formulations  for  AIM 

The  AIM  has  been  used  to  reduce  the  memory  requirement  and  to  accelerate  the  matrix-vector 
multiplication  in  iterative  solver.  The  basic  idea  of  AIM  is  to  approximate  the  far-zone  interac¬ 
tion  using  the  FFT,  and  then  to  compute  the  near-zone  interaction  directly.  This  matrix-vector 
multiplication  can  be  split  into  two  parts,  i.e. 

ZI  =  znear  I  +  zfarl  (2.27) 

where  Z"Lar is  a  sparse  matrix  that  contains  only  the  nearby  elements  within  a  threshold  distance 
and  Z^ar  represents  the  farzone  interaction  of  the  elements. 

The  matrix-vector  multiplication  implemented  inside  the  AIM  solver  can  be  summarized  suing 
the  following  steps  [41]: 


•  project  the  patch  currents,  either  using  multipole  [23]  expansion  or  far-held  matching[42],  to 
the  surrounding  grids; 

•  compute  the  grid  potentials  using  FFT; 

•  interpolate  the  potentials  back  to  the  patch;  and 
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•  compute  the  correct  near- field  interaction  for  each  of  the  patches. 


To  employ  the  AIM,  the  object  is  enclosed  in  a  rectangular  grid  and  then  recursively  subdivided 
into  small  rectangular  grids.  In  order  to  use  FFT  to  approximate  the  far-zone  interaction,  we  need 
to  transform  the  SWG  and  RWG  basis  functions  into  the  Cartesian  grids.  It  is  noted  that  by  using 
vector  identity  and  divergence  theorem,  the  V  operator  in  (2.24)  can  be  shifted  and  operated  on 
the  testing  function  [39],  [40].  Hence  we  can  denote  jn  to  represent  any  one  of  the  components 
y)}  and  express  all  the  matrix  elements  in  (2.24)  in  the  following  unified  form: 

Zmn=  f  f  'ym(r)g(r,r,)^n(r,)dr'dr  (2.28) 

J  Ot-m  ”  Oin 

The  transformation  function  7 n(r)  can  be  approximated  as  a  linear  combination  of  Dirac  delta 
functions,  i.e. 


(m+1)3 

7 n(r)  ~  7 (r)  =  ^2  Anu^r  ~  r')  (2-29) 

U=1 

where  M  is  the  expansion  order  and  Anu  denotes  the  expansion  coefficients  of  'yn(r)-  Anu  can  be 
determined  using  the  multipole  expansion  or  far-held  approximation  [34] .  The  multipole  expansion 
method  is  based  on  the  criterion  that  the  coefficients  Anu  produce  the  same  multipole  moments  of 
the  original  basis  function 


7 n(r)(x  -  *0 )mi(y  -  yo)m2(z  -  z0)m3dr 


(m+1)3 

^2  ( xnu  -  Xo)mi  (: Vnu  ~  J/o)""2  (znu  -  Zo)™3  Ar 


U=  1 


for  0  <  {7711,771.2,171.3}  <  M 


(2.30) 


where  the  reference  point  ro  =  (xo,  yo ,  ^0)  is  chosen  as  the  center  of  the  basis  function.  The  far-held 
approximation  method  matches  the  far-helds  produced  by  yn(r )  and  7 n(r)  on  a  unit  sphere,  and 
it  thus  obtains  the  expansion  coefficients.  Once  the  transformation  function  has  been  determined, 
the  matrix  elements  can  be  approximated  as 

(m+1)3  (771+1)3 

Z’rnn  —  E  E  Am  v9{Xvi  nu  (2.31) 

V=1  U=  1 

By  using  the  transformation  function,  now  we  are  able  to  compute  the  two  components  in  the 
matrix-vector  multiplication  in  A. 12  with 

Zfar  =  AgAT  (2.32a) 

znear  =  z2fM  -  Zfar  (2.32b) 

where  Znear  is  the  matrix  containing  only  the  direct  interaction  of  neighbor  elements  and  A 
represents  the  basis  transformation  matrix  of  the  elements.  The  matrix  g  is  Toeplitz,  and  this 
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enables  the  use  of  FFT  to  compute  the  3-D  convolution  in  the  2.32a  efficiently.  Hence 
represent  the  matrix-vector  multiplication  as 

zi  =  zneari  +  zfari 

Znear/  +  A9-1|9{g}-9{AT7}} 

where  ^{»}  and  9:^1{»}  stand  for  FFT  and  inverse  FFT  respectively. 


we  can 


(2.33) 


Chapter  3 


Numerical  Result  for  AIM 


3.1  Validation  of  the  VSIE-AIM  Algorithm 

An  example  is  considered  for  validating  the  accuracy  of  the  code  developed.  We  consider  a  coated 
dielectric  sphere  having  a  radius  of  1  nr.  The  core  of  the  sphere  has  a  radius  of  0.9  nr  and  a 
relative  permittivity  of  er\  =  1.4  —  0.3j,  and  the  thickness  of  coating  layer  is  O.lnr  with  a  relative 
permittivity  of  er2  =  1.6  —  0.8j.  The  bistatic  radar  cross  section(RCS)  for  the  VV-  and  HH- 
polarizations  at  750  MHz  are  computed  and  the  results  are  shown  in  Figure  3.1.  The  results  are 
compared  with  Mie  series  solutions  and  good  agreements  have  been  observed  in  the  comparisons. 
So,  in  the  next  step  we  apply  this  algorithm  in  the  scattering  problem  of  SRR  structure. 


3.2  Resonance  Frequency  of  SRRs 

3.2.1  Geometry  Statement  of  SRR  structure 

Here  we  choose  square  shape  SRR  structure  instead  of  round  shape  SRR  structure  since  the  it 
is  easier  to  mesh  and  hence,  easier  for  us  to  get  an  accurate  result  since  when  we  mesh  the 
conducting  surface  into  triangular  patches  without  any  mesh-induced  segment  and  apply  the  SWG 
basis  function  to  the  patches. 

3.2.2  Resonance  Frequency  of  SRRs 

First,  we  consider  one  row  of  SRR  structure  (3  inclusions  placed  side  by  side  along  y-).  The 
geometry  and  dimensions  of  this  kind  of  inclusions  are  shown  in  Figure  3.2.  The  width  of  all  metal 
strips  is  0.25nrm,  the  thickness  of  dielectric  is  0.254nrm  and  the  length  of  each  dielectric  square 
is  3. 3mm.  All  the  rings  and  wires  are  place  at  the  center  of  each  dielectric  square.  The  relative 
permittivity  of  dielectric  is  set  to  be  e  =  1. 

We  examine  the  scattering  cross  section  (SCS)  [9]  of  this  row  of  3  elements  and  plot  it  verses 
frequency  in  Figure  3.5.  It  can  be  seen  that  the  resonant  frequency  is  approximately  15.80  GHz. 
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Figure  3.1:  Bistatic  RCS  of  a  coated  dielectric  sphere  (ai  =  0.9  nr,  er\  =  1.4  —  0.3j;  a-i  =  1  m, 
er2  =  1.6  —  0.8j)  at  750  MHz. 


0.46  mm 


Figure  3.2:  The  front  and  back  view  of  1  SRR  inclusion 
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Figure  3.3:  The  rings  view  of  1  row  of  SRRs  (inclusions) 


Figure  3.4:  The  wire  view  of  1  row  of  SRRs  (inclusions) 


frequency  (GHz) 


Figure  3.5:  SCS  versus  frequency  for  a  row  of  inclusions  (3  SRRs)  at  9  =  (f>  =  90' 


Chapter  3.  Numerical  Result  for  AIM 


19 


3.3  RCS  of  a  SRR  Slab  With  Different  Plane- Waves  as  Illumina¬ 
tor 


Next,  at  this  resonant  frequency,  we  analyze  the  propagation  characteristics  of  electromagnetic 
wave  in  a  LHM  sample  placed  in  the  free  space.  It’s  a  LHM  slab  composed  of  many  rows  of 
inclusions  shown  in  Figure  3.2.  The  3-D  view  of  this  slab  is  depicted  in  3.6.  The  space  distances  of 
the  inclusions  denoted,  respectively,  by  dx,  dy  and  dz  in  the  x—,  y—  and  z—  directions,  are  all  3.3 
mm.  There  are  totally  54  SRR  elements  which  are  arranged  into  18  rows  along  x—  and  3  columns 
along  y—.  In  order  to  show  that  our  AIM  algorithm  is  suitable  for  this  structure,  we  illuminate 
the  slab  with  two  different  plane  waves  coming  from  —  y  direction  as  shown  in  Figure  3.6,  and  then 
we  check  the  RCS  of  this  SRR  slab  in  each  case.  Actually,  SRR  requires  the  electric  field  Einc  in 
parallel  with  the  plane  of  the  ring  ( yOz  plane)  to  get  maximum  magnetic  resonant.  In  other  words, 
the  rings  are  not  supposed  to  be  on  the  H  —  k  plane  [14] .  So  if  Ernc  is  perpendicular  with  the  plane 
of  rings  as  show  in  case  <  2  >,  the  RCS  must  be  very  small  and  the  whole  structure  becomes  no 
more  than  an  ordinary  scatterer  which  does  not  acquire  negative  refractive  index.  Figure  3.3  to 
Figure  3.10  show  values  of  Ez  and  RCS  when  the  illuminators  are  of  case  <  1  >  and  case  <  2  > 
respectively. 


Figure  3.6:  Structure  of  the  slab  composed  by  SRRs  and  wires  with  and  two  different  incidents  in 
the  reference  frame 

Another  advantage  of  our  program  compared  with  other  algorithm  is  that  we  include  volume 
element  here.  Thus,  we  can  change  the  e  to  observe  the  change  of  Ez  and  RCS  very  conveniently. 
Here  we  use  e  =  1.60  —  0.05j  instead  of  e  =  1.  Below  are  the  plot  of  Ez  and  RCS  versus  </>  at  15.8 
GHz  when  the  input  are  of  case  <  1  >  and  case  <  2  >  respectively. 

Further  more,  Figures  3.15-3.18  listed  below  show  the  Ez  and  RCS  of  e  =  2.40  —  0.20j  at 
15.8  GHz  when  the  slab  under  magnetic  resonance  and  no  magnetic  resonance  input  respectively. 
Before  calculating  the  scattering  field  of  this  material,  we  first  examine  the  resonant  frequency 
again  to  ensure  the  new  slab  works  in  the  pass  band  of  resonance.  We  found  that  the  central 
frequency  changes  less  than  0.1  Ghz  which  is  actually  a  very  small  number  compared  with  the 
central  frequency.  So,  we  can  still  fix  the  frequency  at  15.8  GHz  which  is  in  the  pass  band.  From 
the  Ez  and  RCS  result,  we  can  observe  a  decrease  in  the  scattering  field  which  is  due  to  higher 
permittivity. 
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Figure  3.7:  Ez  of  the  slab  of  case  <  1  >  at  15.8  GHz 


o 


Figure  3.8:  RCS  of  the  slab  of  case  <  1  >  at  15.8  GHz 
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Figure  3.9:  Ez  of  the  slab  of  case  <  2  >  at  15.8  GHz 
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Figure  3.10:  RCS  of  the  slab  of  case  <  2  >  at  15.8  GHz 
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Figure  3.11:  Ez  of  the  slab  of  case  <  1  >  at  15.8  GHz  when  e  =  1.60  —  0.05j 
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Figure  3.12:  RCS  of  the  slab  of  case  <  1  >  at  15.8  GHz  when  e  =  1.60  —  0.05j 
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Figure  3.13:  Ez  of  the  slab  of  case  <  2  >  at  15.8  GHz  when  e  =  1.60  —  0.05j 
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Figure  3.14:  RCS  of  the  slab  of  case  <  2  >  at  15.8  GHz  when  e  =  1.60  —  0.05j 
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Figure  3.15:  Ez  of  the  slab  of  case  <  1  >  at  15.8  GHz  when  e  =  2.40  —  0.20 j 


Figure  3.16:  RCS  of  the  slab  of  case  <  1  >  at  15.8  GHz  when  e  =  2.40  —  0.20 j 
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Figure  3.17:  Ez  of  the  slab  of  case  <  2  >  at  15.8  GHz  when  e  =  2.40  —  0.20 j 
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Figure  3.18:  RCS  of  the  slab  of  case  <  2  >  at  15.8  GHz  when  e  =  2.40  —  0.20 j 
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3.4  Conclusion  for  AIM  Algorithm 

The  structures  of  SRR  with  wires  behave  as  a  magnetic  conductor  in  the  vicinity  of  resonant 
frequency.  By  applying  the  fast  solver  AIM,  the  CPU  time  per  iteration  as  well  as  memory 
requirement  is  greatly  reduced.  Even  though  some  very  expensive  EM  softwares  can  calculate 
RCS  and  Ez  for  small  scale  problems,  they  are  not  capable  in  simulating  large  problems  such  as 
SRR  slab.  SRR  structures  are  often  electrically  small  but  have  a  large  number  of  unknowns,  so 
AIM  is  especially  suitable  to  analyze  such  kind  of  LHMs.  And,  since  SRR  is  composed  of  dielectric 
body  and  metal  surface,  we  introduce  both  volume  and  surface  into  our  program.  After  a  rough 
validation  of  the  code,  we  first  use  a  few  inclusions  to  find  the  resonant  frequency,  then  at  the 
vicinity  of  this  frequency,  an  SRR  slab  is  tested  when  the  illuminator  are  two  different  plane  waves 
while  only  the  first  one  can  induce  magnetic  resonance.  We  find  that  when  the  magnetic  resonance 
is  achieved,  most  power  can  be  transferred  to  the  other  side  of  the  slab  which  satisfies  the  property 
of  LHM.  However,  when  there  is  no  magnetic  resonance,  most  power  will  not  be  able  to  go  through 
the  slab.  Then  we  change  the  permittivity  of  dielectric  volume  to  observe  it’s  effect  on  scattered 
field.  The  result  prove  that  we  achieve  high  precision  with  the  lowest  memory  and  CPU  time  cost, 
and  that  our  algorithm  is  efficient  enough  for  any  composite  LHM  design  and  simulation. 
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Application  of  Left-handed  Material 
with  New  Designs 

4.1  Microstrip  Lines  Loaded  With  Flower-shaped  Left-handed  Res¬ 
onators 

4.1.1  Introduction  of  Microstrip  Lines  Loaded  with  CSRR  and  Capatitors 

In  this  work,  microstrip  lines  loaded  with  flower-shaped  left-handed  resonators  and  series  inter¬ 
digital  capacitors  are  modeled  by  means  of  a  lumped  element  equivalent  circuit.  The  structures 
exhibit  a  left-handed  (or  backward  wave)  behavior  at  a  relatively  broad  frequency  range.  That  is  a 
huge  improvement  to  the  conventional  SRR  structure  which  has  very  narrow  bandwidth  and  high 
loss.  Also,  we  explored  the  equivalent  circuit  at  the  left-handed  frequency.  The  equivalent  circuit 
model  is  able  to  accurately  describe  the  frequency  responses  of  the  structures  in  the  left-handed 
transmission  band,  the  forbidden  band,  and  the  lower  frequency  region  of  the  forward  allowed 
band.  The  limitations  of  this  circuit  model  to  properly  predict  the  behavior  of  the  structures  at 
higher  frequencies  (mainly  related  to  the  Bragg  effect  caused  by  periodicity)  are  pointed  out. 

Complementary  split  rings  resonators  (CSRRs)  have  been  recently  proposed  by  some  of  the 
authors  of  this  manuscript  as  new  resonant  constitutive  elements  for  the  synthesis  of  negative 
permittivity  nretamaterials  [43].  By  combining  these  resonators  with  series  capacitive  gaps,  one- 
dimensional  left-handed  structures  were  implemented  in  microstrip  technology  [44]. The  required 
negative  value  of  the  effective  permeability  to  achieve  left-handedness  was  due  to  the  presence  of 
the  series  gaps,  which  make  the  structure  to  behave  as  a  magnetic  plasma,  with  negative  value  of 
permeability  up  to  a  cutoff  (or  plasma)  frequency.  The  frequency  response  of  CSRR/gap-loaded 
nricrostrip  lines  exhibits  a  left-handed  (or  backward  wave)  transmission  band  in  that  region  were 
negative  permeability  and  permittivity  simultaneously  coexist,  and  this  typically  occurs  in  a  narrow 
(or  moderate)  band  below  the  plasma  frequency  of  the  gap-related  structure,  i.e.  in  the  vicinity 
of  the  resonance  frequency  of  CSRRs  [44].  Above  this  band,  the  effective  permittivity  changes 
sign  and  signal  propagation  is  inhibited.  However,  above  the  plasma  frequency  also  the  effective 
permeability  of  the  artificial  line  switches  to  a  positive  value,  and  signal  propagation  is  allowed, 
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although  wave  transmission  is  forward,  rather  than  backward.  Thus,  the  typical  behavior  of  these 
resonant  type  left-handed  lines  is  indeed  similar  to  that  of  the  so-called  composite  right- /left- 
handed  (CRLH)  transmission  lines,  where  a  left-handed  and  a  right-handed  transmission  bands, 
separated  by  a  frequency  gap,  are  present  [19].  In  CRLH  lines,  the  left-handed  band,  at  lower 
frequencies,  is  due  to  the  presence  of  series  gaps  and  shunt  stubs  loading  the  line,  whereas  forward 
wave  transmission  is  related  to  the  parasitic  elements  of  the  host  line.  To  properly  model  the 
behavior  of  CSRR-based  left-handed  microstrip  lines,  a  lumped  element  equivalent  circuit  model 
has  been  previously  proposed  [45].  This  model  has  been  found  to  describe  the  frequency  response 
in  the  left-handed  band  to  a  good  approximation.  And  the  capability  of  the  model  to  also  predict 
transmission  line  behavior  above  that  band  is  explored  by  [46]. 


4.1.2  The  Novel  Planar  Negative  Refractive  Index  Structure 

A  novel  2D  planar  distributed  NRI  structure  which  is  very  likely  the  performance  of  CSRR,  is 
proposed  previously  by  Naoko  Natrunaga  etc  [47].  The  structure  consists  of  a  2D  periodic  array 
of  unit  cells  with  metal  patterns  on  the  both  sides  of  a  substrate  and  it  does  not  require  any  vias. 
Without  vias,  is  the  structure  scalable  to  millimeter-wave  frequency  range  or  above.  It  is  also 
advantageous  that  the  structure  can  be  fabricated  by  using  conventional  MMIC  processes.  The 
structure  were  experimentally  verified  of  the  left-handed  characteristics. 

The  structure  is  the  planar  structure  with  metal  patterns  on  both  the  upper  and  lower  layers 
of  the  substrate  as  shown  in  Figure  4.1(a)  and  (b),  respectively.  Let  us  consider  a  square  unit  cell 
with  the  period  a  in  both  the  x  and  y  directions  for  simplicity.  The  unit  cell  is  depicted  with  the 
dotted  line  in  Figure  4.1(a)  and  (b).  On  the  upper  layer,  are  the  45-  degree  rotated  square  metal 
patches  with  the  sides.  On  the  lower  layer,  four  isosceles  triangles  surrounded  by  a  square  frame 
are  connected  at  their  apex,  and  from  the  connected  portion,  four  narrow  metal  strips  are  radially 
connected  to  the  outer  frame  in  the  unit  cell.  The  upper  layer  pattern  lies  on  the  lower  layer 
pattern  with  the  relation  shown  in  4.1(c).  The  upper  metal  patch  overlaps  the  isosceles  triangles 
forming  an  MIM  capacitor,  which  provides  a  large  series  capacitance  between  the  adjacent  unit 
cells  and  a  frequency  band  of  negative  e.  The  four  radial  strips  bridge  the  series  capacitor  and  the 
outer  frame.  When  the  area  of  the  outer  frame  is  large  enough  and  has  a  quasi-ground  voltage,  the 
bridging  strips  work  as  shunt  inductances  and  can  provide  negative  p  at  a  frequency  band.  When 
the  series  capacitance  and  the  shunt  inductance  are  large  enough,  in  other  words,  the  negative  e 
and  negative  p  appear  at  the  same  frequency  band,  the  structure  is  verified  that  it  can  work  as  a 
left-handed  structure  during  a  particular  frequency  [47] . 

However,  before  utilizing  the  structure  in  microstrip  line,  we  find  the  bandwidth  of  such  structure 
as  shown  in  Fig  4.2can  be  further  improved  by  the  changing  the  shape  of  isosceles  triangles  into 
the  flower-shape  strip  as  shown  in  4.3.  To  prove  with  this  shape  can  obtain  better  performance,  we 
etch  these  two  structures  into  the  ground  plane  of  a  microstrip  line  to  see  the  transmission  ability 
vs  frequency.  As  we  stated  before,  the  isosceles  triangle  can  provide  negative  p  in  a  frequency  band 
which  will  be  stopband  if  we  want  to  observe  the  S'21  of  these  microstrip  lines. 

After  we  examined  the  stop-band  of  the  structures  shown  in  Figure  4.2  and  Figure  4.3.  The 
stop-band  for  the  original  triangle  structure  is  from  7-8GHz  as  shown  in  Figure  4.4,  while  the 
modified  flower  structure’s  stop-band  is  from  8.8-10.6GHz  as  shown  in  Figure  4.5.  So  we  add 
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Figure  4.1:  Proposed  2D  NRI  structure,  (a)  Upper  layer,  (b)  Lower  layer,  (c)  Unit  cell. 


Figure  4.2:  Top  view  of  a  microstrip  line  with  2  rows  of  triangle-shaped  left-handed  pattern  on  the 
ground  plane 
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Figure  4.3:  Top  view  of  a  microstrip  line  with  2  rows  of  flower-shaped  left-handed  pattern  on  the 
ground  plane 
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Figure  4.4:  S-parameter  of  a  microstrip  line  with  2  row  of  triangle-shape  left-handed  pattern  on 
the  ground  plane  and  interdigital  capacitors  on  the  stripline  (red  line  represents  S21  and  blue  line 
represents  <Sii 

capacitor  onto  the  microstrip  line(Figure  4.6  and  Figure  4.7),  which  finally  provide  a  C-L  circuit 
and  induce  both  negative  p  and  negative  e.  Now  the  stop-band  becomes  a  pass-band,  and  most 
important,  the  pass-band  is  as  wide  as  2.5GHz  as  shown  in  Figure  4.8. 


4.2  Microstrip-feed  Patch  Antenna  with  Left-handed  Elements 

4.2.1  Novel  Patch  Antenna  with  LH  Element  and  Its  Performance 

In  this  section,  we  apply  the  triangle  structure  into  the  fabrication  of  patch  antenna.  The  model 
we  designed  is  shown  in  Figure  4.9  and  Figure  4.10.  And  the  Si  \  of  the  proposed  antenna  is  shown 
in  Figure  4.11.  We  can  see  that  the  bandwidth  is  from  5.1-9.2GHz.  The  VSWR  is  shown  in  Figure 
4.12.  Obviously,  during  the  working  band,  the  VSWR  is  effectively  suppressed  below  2.  The  3D 
radiation  pattern  at  7GHz  and  8.1GHz  is  depicted  in  Figure  4.14  and  Figure  4.15  respectively. 
We  can  see  that  the  radiation  is  mainly  gathering  at  the  direction  of  8  =  90° ,  cj)  =  90°  which  is 
quite  different  from  a  traditional  patch  antenna’s  main  radiation  direction  (9  =  0°).  This  can  be 
explained  by  the  transmission  of  backward  wave  from  the  feed  line  to  the  end  of  the  substrate. 

At  the  frequency  of  7GHz,  on  9  =  90°  plane,  as  shown  in  Figure  4.16,  the  main  lobe  magnitude 
of  the  antenna  is  4.9dBi.  The  main  lobe  direction  is  at  <f>  =  105°  and  the  angular  width  (3dB)  is 
67.1°.  The  side  lobe  level  is  — 3.5dB. 

At  the  frequency  of  8.1GHz,  on  9  =  90°  plane,  as  shown  in  Figure  4.17,  the  main  lobe  magnitude 
of  the  antenna  is  6.3dBi.  The  main  lobe  direction  is  at  </>  =  100°  and  the  angular  width  (3dB)  is 
57.2°.  The  side  lobe  level  is  — 4.8dB. 

At  the  frequency  of  7GHz,  on  <j>  =  90°  plane,  as  shown  in  Figure  4.18,  the  main  lobe  magnitude 
of  the  antenna  is  4.2dBi,  which  is  quite  favorable  for  a  single  patch  antenna,  the  main  lobe  direction 
is  at  (j)  =  80°  and  the  angular  width  (3dB)  is  142.8°. 
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Figure  4.5:  S-parameter  of  a  microstrip  line  with  2  row  of  flower-shape  left-handed  pattern  etched 
on  the  ground  plane  and  interdigital  capacitors  on  the  stripline  (red  line  represents  S21  and  blue 
line  represents  <Sii 


Figure  4.6:  Top  view  of  a  microstrip  line  with  1  row  of  flower-shape  left-handed  pattern  on  the 
ground  plane  and  interdigital  capacitors  on  the  stripline 


Figure  4.7:  Bottom  view  of  a  microstrip  line  with  1  row  of  flower-shape  left-handed  pattern  on  the 
ground  plane  and  interdigital  capacitors  on  the  stripline 
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Figure  4.8:  S-parameter  of  a  microstrip  line  with  1  row  of  flower-shape  left-handed  pattern  etched 
on  the  ground  plane  and  interdigital  capacitors  on  the  stripline  (red  line  represents  S21  and  blue 
line  represents  <Sii 


Figure  4.9:  Top  view  of  the  simulated  antenna  (blue  part  is  substrate  while  pink  part  shows  metal 
surface) 
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Figure  4.10:  Bottom  view  of  the  simulated  antenna  (blue  part  is  substrate  while  pink  part  shows 
metal  surface) 
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Figure  4.11:  Simulated  Sn  against  frequency  of  patch  antenna  with  triangle  shaped  gaps  on  patch 
and  strip  bridges  on  ground 
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Figure  4.12:  Simulated  VSWR  against  frequency  of  patch  antenna  with  triangle  shaped  gaps  on 
patch  and  strip  bridges  on  ground 
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Figure  4.13:  Simulated  antenna  balance  vs  frequency  in  the  working  band 
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Figure  4.14:  Simulated  3D  radiation  pattern  at  the  frequency  of  7GHz 
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Figure  4.15:  Simulated  3D  radiation  pattern  at  the  frequency  of  8.1GHz 
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Figure  4.16:  Simulated  radiation  pattern  at  7GHz  on  the  plane  of  6  =  90° 
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Figure  4.17:  Simulated  radiation  pattern  at  8.1GHz  on  the  plane  of  6  =  90' 
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Figure  4.18:  Simulated  radiation  pattern  at  7GHz  on  the  plane  of  cj)  =  90° 


At  the  frequency  of  8.1GHz,  on  <f>  =  90°  plane,  as  shown  in  Figure  4.19,  the  main  lobe  magnitude 
of  the  antenna  is  5.8dBi.  The  main  lobe  direction  is  at  6  =  90°  and  the  angular  width  (3dB)  is 
114.4°.  The  side  lobe  level  is  about  — 7.6dB. 


4.2.2  Conclusion  for  the  Application 

We  fabricated  the  patch  antenna  as  shown  in  Figure  4.20  and  Figure  4.21  with  left-handed  patterns 
and  gaps  etched  on  the  patch  and  ground  plane  respectively.  We  find  a  good  agreement  between 
the  experimental  results  and  the  simulation  prediction.  This  antenna  overcome  the  disadvantage  of 
traditional  patch  antenna’s  narrow  bandwidth  by  enhance  the  radiation  with  left-handed  element 
we  present  in  the  previous  sections.  The  radiation  efficiency  is  very  high,  which  shows  most  energy 
is  radiated  rather  than  wasted  thermally.  At  the  mean  time,  the  gain  of  such  LH-enhanced  antenna 
is  high,  and  the  side  lobes  are  well  suppressed. 
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Figure  4.19:  Simulated  radiation  pattern  at  8.1GHz  on  the  plane  of  <j)  =  90° 


Figure  4.20:  Top  view  of  the  fabricated  antenna  (£=2.2,  substrate  thickness  h=0. 794mm,  substrate 
size=28.1mmx32mm,  patch  size=12mmxl6mm) 
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Figure  4.21:  Bottom  view  of  the  fabricated  antenna  (width  of  gaps=0.4mm) 
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Base  on  previous  result,  we  can  see  an  obvious  pass  band  of  the  SRR  slab.  To  be  more  convincing, 
we  are  going  to  show  more  LH  property  of  this  SRR  structure  such  as  the  negative  refractive  index. 

In  our  future  work,  we  plan  to  build  a  SRR  prism  as  shown  in  Figure  5.2  and  Figure  5.3.  We  can 
use  the  devices  shown  in  Figure  5.1  to  measure  the  angle  of  refraction.  The  prism  we  design  still 
uses  the  basic  SRR  element  shown  in  Figure  3.2.  The  whole  prism  is  shown  in  Figure  5.2  and  Figure 
5.3.  This  model  is  built  and  meshed  in  FEKO  5.0  which  require  9G  memory  to  apply  traditional 
MoM.  However,  with  our  AIM  algorithm,  it  is  estimated  that  only  2G  memory  is  required. 

With  this  structure,  we  are  expecting  to  observe  an  obvious  negative  refraction  from  the  prism. 
If  this  simulation  can  be  performed  successfully,  we  will  go  on  with  the  implementation  of  our  code. 
For  example,  we  will  adopt  periodic  Green’s  function  to  further  reduce  computation  complexity.  As 
for  the  design  of  antenna,  we  can  also  include  multi-layer  Green’s  function  to  enhance  the  accuracy 
as  well  as  the  computation  speed  of  the  code  [48]  [49].  Because  a  full-wave  Greens  function  com¬ 
pletely  define  the  field  inside  a  multi-layer  dielectric  structure  due  to  a  current  element  arbitrarily 
placed  between  any  two  layers  in  two-dimensional  spectral-domain  form.  It  is  derived  by  solving  a 
’’standard”  form  containing  the  current  element  with  two  substrates  on  either  side  of  it,  and  using 
an  iterative  algorithm  to  take  care  of  additional  layers.  Another  iterative  algorithm  is  then  used  to 
find  the  field  in  any  layer  in  terms  of  the  field  expressions  in  the  two  layers  of  the  ’’standard”  form. 


Microwave  absorber 


Figure  5.1:  Special  device  for  measurements  of  angles  of  refractions 
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Figure  5.2:  3D  view  of  SRR  prism 


Figure  5.3:  Top  view  of  SRR  prism 
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The  locations  of  the  poles  of  the  Greens  function  can  be  predicted.  Higher  order  interpolatory  basis 
functions  defined  on  curvilinear  triangular  patches  also  can  be  used  to  provide  necessary  flexibility 
and  accuracy  for  the  discretization  of  arbitrary  shapes  and  to  offer  a  better  convergence  than  lower 
order  basis  functions.  The  combination  of  the  improved  DCIM  and  the  higher  order  basis  functions 
results  in  an  efficient  and  accurate  MoM  analysis  for  3-D  multi-layer  microstrip  structures.  So  this 
kind  of  Green’s  function  can  be  applied  to  a  few  multilayer  transmission  or  radiation  structures. 
After  all  these  work,  our  algorithm  will  become  an  even  more  powerful  numerical  simulation  tool. 

As  to  the  design  of  new  left-handed  structure  and  application  aspect,  we  can  first  extent  the 
LH-enhanced  patch  antenna  into  antenna  array  and  see  what  further  improvement  can  be  made. 
Then  we  can  adopt  various  structures  discovered  with  LH  characteristics  into  a  large  variety  of 
devices  such  as  microstrip  filter,  power  divider,  coupler  [50]  as  well  as  antenna  fabrication  [51]. 
Many  theoretical  predictions  and  new  LH  patterns  design  have  been  done  to  serve  the  application. 
However,  to  date,  the  design  and  fabrication  of  practical  devices  is  still  mainly  based  on  the 
traditional  SRR  or  H  structure  which  has  extremely  high  loss  and  narrow  bandwidth.  To  open  a 
new  world  of  LHM  application,  we  are  going  to  make  better  use  of  current  research  result  and  try 
to  explore  various  practical  applications. 
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A.l  RWG  Basis  Function  —  Triangular  Modeling  Method  for  Sur¬ 
face  Element 

A. 1.1  Development  of  RWG  Basis  Function 

Let  S  denote  the  surface  of  an  open  or  closed  perfectly  conducting  scatterer.  SWG  is  the  short 
for  Rao-Wilton-Glisson  basis  function.  Each  basis  function  is  associated  with  a  interior  edge  (i.e., 
nonboundary  edge  )  of  the  patch  model  and  is  to  vanish  everywhere  on  S  except  in  the  two  triangles 
attached  to  that  edge.  Figure  A.l  shows  two  such  triangles,  T+  and  T~,  corresponding  to  the  nth 
edge  of  triangulated  surface  modeling  a  scatterer.  Points  in  T+  may  be  designated  either  by  the 
position  vector  r  defined  with  respect  to  O,  or  by  the  free  vertex  of  T+.  Similar  remarks  apply  to 
the  position  vector  p~  except  that  it  is  directed  toward  the  free  vertex  of  T~ .  The  plus  or  minus 
designation  of  the  triangles  is  determined  by  the  choice  of  a  positive  current  reference  direction  for 
the  nth  edge,  the  reference  for  which  is  assumed  to  be  from  T, ’+  to  T~ .  We  define  the  vector  basis 
function  associated  with  the  nth  edge  as 


Figure  A.l: 


Triangle  pair  and  geometrical  parameters  associated  with  interior  edge 
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Figure  A. 2:  Geometry  for  construction  of  component  of  basis  function  normal  to  edge 
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(A.l) 


where  ln  is  the  length  of  the  edge  and  is  the  area  of  triangle  T±.  (Note  that  we  use  the 
convention,  followed  throughout  the  report,  that  subscripts  refer  to  edges  while  superscripts  refer 
to  faces.)  The  basis  function  fn  is  used  to  approximately  represent  the  surface  current,  and  we 
list  and  discuss  below  some  properties  which  make  it  uniquely  suited  to  this  role. 

1)  The  current  has  no  component  normal  to  the  boundary  (which  excludes  the  common  edge) 
of  the  surface  formed  by  the  triangle  pair  T+  and  T“,  and  hence  no  line  charges  exist  along  this 
boundary. 

2)  The  component  of  current  normal  to  the  nth  edge  is  constant  and  continuous  across  the 
edge  as  may  be  seen  with  the  aid  of  Figure  A. 2,  which  shows  that  the  normal  component  of  p± 
along  edge  n  is  just  the  height  of  triangle  with  edge  n  as  the  base  and  the  height  expressed  as 
(2 A±)/ln.  This  latter  factor  normalizes  fn  in  such  that  its  flux  density  normal  to  edge  n  is  unity, 
ensuring  continuity  of  current  normal  to  the  edge.  This  result,  together  with  1),  implies  that  all 
edges  of  T+  and  T~  are  free  of  line  charges. 

3)  The  surface  divergence  of  fn ,  which  is  proportional  to  the  surface  charge  density  associated 
with  the  basis  element,  is 


Vs  •  /n(r) 


0, 
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(A.2) 


since  the  surface  divergence  in  is  ^  •  The  charge  density  is  thus  constant  in  each 

triangle,  the  total  charge  associated  with  the  triangle  pair  T+  and  T~  is  zero,  and  the  basis 


functions  for  the  charge  evidently  have  the  form  of  pulse  doublets. 
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4)  The  moment  of  fn  is  given  by  (A+  +  An)f™9  where 

(A+  +  A^)fr=  fT.n++T-fndS 

=  lf(Pn++Pn~) 

=  ln{rc+  +  rc~)  (A. 3) 

and  pc±  is  the  vector  between  the  free  vertex  and  the  centroid  of  with  p^~  directed  toward 
and  p^f  directed  away  from  the  vertex,  as  shown  in  Figure  A.l,  and  is  the  vector  from  O  to 
the  centroid  of  T±.  (A. 3)  may  be  most  easily  derived  by  expressing  the  integral  in  terms  of  area 
coordinates,  to  be  discussed  below. 

The  current  on  S  may  be  approximated  in  terms  of  the  fn,  as 

N 

J^^Infn(r)  (A.4) 

71=  1 

where  N  is  the  number  of  interior  (nonboundary)  edges.  Since  a  basis  function  is  associated 
with  each  nonboundary  edge  of  the  triangulated  structure,  up  to  three  basis  functions  may  have 
nonzero  values  within  each  triangular  face.  But  at  a  given  edge  only  the  basis  function  associated 
with  that  edge  has  a  current  component  normal  to  the  edge  since,  according  to  1),  all  other  basis 
currents  in  adjacent  faces  are  parallel  to  the  edge.  Furthermore,  since  the  normal  component 
Of  fn,  at  the  nth  edge  is  unity,  each  coefficient  In,  in  (A.4)  may  be  interpreted  as  the  normal 
component  of  current  density  flowing  past  the  nth  edge.  Also,  we  see  that  the  basis  functions  are 
independent  in  each  triangle  since  the  current  normal  to  the  nth  edge,  In  in  (A.4),  is  an  independent 
quantity.  At  surface  boundary  edges,  the  sum  of  the  normal  components  of  current  on  opposite 
sides  of  the  surface  cancel  because  of  current  continuity.  Therefore  we  neither  define  nor  include  in 
(A.4)  contributions  from  basis  functions  associated  with  such  edges.  Because  of  the  considerable 
variation  in  the  direction  of  the  flow  lines  of  f  n  within  a  triangle,  it  is  not  at  first  obvious  that  a 
linear  superposition  of  basis  functions  is  capable  of  representing,  say,  a  constant  current  flowing 
in  an  arbitrary  direction  within  a  triangle.  That  this  is  possible,  however,  can  be  seen  with  the 
aid  of  Figure  A. 3,  which  shows  a  triangle  Tq  with  edges  arbitrarily  labeled  1,  2,  and  3  (in  effect, 
we  employ  here  a  local  indexing  scheme,  in  contrast  to  the  global  indexing  scheme  used  earlier). 
With  the  vectors  pi,  p2  and  pz  as  shown,  the  basis  functions  in  Tq  are  =  ( li/2Aq)p.i ,  i  =  1, 
2,  3,  where  Aq  is  the  triangle  area  and  where,  for  simplicity,  the  current  reference  directions  are 
assumed  to  be  out  of  the  triangle  for  each  edge.  It  is  apparent  from  the  figure  and  the  definition 
of  fi  that  the  linear  combinations  Z2/1  —  ffr2  and  hf  1  ~  hfs  are  constant  vectors  for  every  point 
r  in  Tq  and  are  parallel  to  sides  3  and  2,  respectively.  Since  the  two  composite  forms  are  linearly 
independent  (i.e. ,  nonparallel),  a  constant  vector  of  arbitrary  magnitude  and  direction  within  Tq 
may  be  synthesized  by  an  appropriate  linear  combination  of  the  two  forms,  as  asserted. 


A. 1.2  Testing  Procedure  of  RWG  Basis  Function 

The  next  step  in  the  method  of  moments  is  to  select  a  testing  procedure.  We  choose  as  testing 
functions  the  expansion  functions  fn,  developed  in  the  previous  section.  With  a  symmetric  product 
defined  as 
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Figure  A. 3:  Local  coordinates  and  edges  for  source  triangle  Tq  with  the  observation  point  in 
triangle  Tp 


<  f,g>=  [  f  ■  gds  (A. 5) 

Js 

we  derive  an  integrodifferential  equation  for  J  by  enforcing  the  boundary  condition  h  x  ( E{nc  + 
Esca)  =  0  on  S',  obtaining 

-EtZ  =  iy—jioA  -  V$)tan,  r  on  S  (A.6) 

Equation  (A.6)  is  tested  with  fm ,  yielding 

<  Einc,  fm  >=  ju  <  A,  fm  >  +  <  V4>,  fm  >  (A.7) 

If  one  makes  use  of  a  surface  vector  calculus  identity  and  the  properties  of  f  m  at  the  edges  of  S, 
the  last  term  in  (A.7)  can  be  rewritten  as: 


<  V$,/m>=-  [  4>VS  •  fmdS 
Js 

With  (A. 2),  the  integral  in  (A. 8)  may  now  be  written  and  approximated  as  follows: 


$VS  •  fmdS 


=  l 
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A  m  Jt+  An 

=  lmMrc+)  -  3>(rc-)] 


(A.8) 


(A.9) 


In  (A.9)  the  average  of  over  each  triangle  is  approximated  by  the  value  of  $  at  the  triangle 
centroid.  With  similar  approximations,  the  vector  potential  and  incident  field  terms  in  (A.7)  may 
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be  written  as 
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where  the  integral  over  each  triangle  is  eliminated  by  approximating  Einc  (or  A)  in  each  triangle 
by  its  value  at  the  triangle  centroid  and  carrying  out  integrations  similar  to  those  used  to  obtain 
(A. 3).  With  (A. 2  A. 9  A. 10),  (A.7)  now  becomes: 


jiolm 

— 


A(rc+)  ■  &  +  A(rcm)  ■  ?f 
Einc(rc+)-?f  +  Einc(rc-)-^- 


(A- 11) 


which  is  the  equation  enforced  at  each  triangle  edge,  m  =  1,  2,  ...IV. 


We  remark  that  another  interpretation  of  the  testing  procedure  arriving  at  A.  11  is  also  possible. 
One  may  equate  line  integrals  of  the  form  jc  F  -dr,  where  F  represents  the  right  and  left  sides  of 
(A. 6),  and  Crn  is  the  piecewise  linear  path  from  the  point  r to  the  midpoint  of  edge  nr  and  thence 
to  r ^  ■  Einc  and  A  can  be  approximated  along  each  portion  of  the  path  by  their  respective  values 
at  the  triangle  centroids.  The  resulting  equality,  apart  from  the  factor  lm ,  is  (A. 11).  Under  either 
interpretation,  the  testing  procedure  reduces  the  differentiability  requirement  on  $  in  (A. 6)  by 
integrating  V4>  first,  the  procedure  having  been  constructed  with  this  goal  in  mind.  The  purpose 
of  approximations  (A. 9)  and  (A. 10)  is  to  eliminate  surface  integrals  of  the  potential  quantities, 
allowing  a  double  surface  integral  to  be  approximated  by  a  quantity  involving  a  single  surface 
integral  in  the  numerical  computation  of  the  moment  matrix  elements.  These  approximations  are 
justified  by  observing  that  the  potentials  are  locally  smooth  within  each  subdomain,  as  follows 
from  their  integral  definitions  and  the  locally  smooth  nature  of  the  source  representation  in  terms 
of  the  basis  functions. 


A. 1.3  Elements  of  the  Matrix  and  Excitation  Vector  of  RWG 

Substitution  of  the  current  expansion  (A. 4)  into  (A.  11)  yields  an  N  x  N  system  of  linear  equations 
which  may  be  written  in  matrix  form  as 


ZI  =  V 


(A.  12) 


where  Z  =  [Zmn]  is  an  N  x  N  matrix  and  /  =  [In]  and  V  =  [Um]  are  column  vectors  of  length  N. 
Elements  of  Z  and  V  are  given  by 


Zrnin  —  lr 
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Pm 
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1  mn  mn 


(A.13) 
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where 


Pm  i  in—  Pm 


Vm  =  lm  E^-^  +  E 
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4vr  Js  Rm 


o  jkRn 
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^  mn 
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R±  =  I  rc±  -  r'\ 

m  r  m  '  I 
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± 
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Einc(r 


c±) 

m  ) 


(A.14) 

(A.15) 

(A.16) 

(A.17) 

(A.18) 


For  plane  wave  incidence,  we  set 

Ei(r)  =  ( Eee0  +  Effr 0)ejkr  (A.19) 

where  the  propagation  vector  k  is 

k  =  k( sin  9q  cos  fox  +  sin  6q  sin  foy  +  cos  9qz)  (A.20) 

and  (00)  fo)  defines  the  angle  of  arrival  of  the  plane  wave  in  terms  of  the  usual  spherical  coordinate 
convention.  Unit  vectors  $o  and  fo  are  constant  vectors  which  coincide  with  the  usual  spherical 
coordinate  unit  vectors  only  at  points  on  the  line  from  O  in  the  direction  of  k. 

Once  the  elements  of  the  moment  matrix  and  the  forcing  vector  V  are  determined,  one  may  solve 
the  resulting  system  of  linear  equations  (A.  12)  for  the  unknown  column  vector  I.  The  elements  of 
Z  in  (A. 12)  may  be  evaluated  by  naively  computing  Zmn  directly  by  (A. 13)  (with  aid  of  (A.15) 
and  (A.16))  for  each  index  combination  m  and  n.  However,  as  shown  in  the  next  section,  this 
procedure  is  extremely  inefficient  since  the  integrals  required  for  each  combination  of  m  and  n  are 
generally  required  for  a  number  of  other  combinations  as  well. 


A. 2  SWG  Basis  Function  —  Tetrahedral  Modeling  Method  for 
Volume  Element 

A. 2.1  Development  of  SWG  Basis  Function 

The  volume  V  is  assumed  to  be  subdivided  into  a  number  of  tetrahedral  elements  such  that  a 
homogeneous  or  inhomogeneous  dielectric  region  is  approximated  by  a  number  of  tetrahedrons,  in 
each  of  which  the  dielectric  properties  are  approximated  as  constant.  A  homogeneous  dielectric 
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Figure  A. 4:  Pair  of  tetrahedrons  and  geometrical  parameters  associated  with  nth  face 


region  is  bounded  by  a  surface  that  is  approximated  by  triangular  faces  assigned  so  as  to  fit  the 
shape  of  the  surface  well. 


Once  the  volume  of  the  scatterer  has  been  appropriately  modeled  by  tetrahedral  volumes,  the 
faces  of  the  tetrahedrons  are  of  primary  importance  for  the  development  of  the  basis  functions. 
Figure  A.4  shows  two  tetrahedrons  T and  T~,  associated  with  the  nth  face  of  the  subdivided 
region  V  modeling  a  scatterer.  Points  in  T+  may  be  designated  either  by  the  position  vector 
position  vector  r  defined  with  respect  to  O,  or  by  the  position  vector  defined  with  respect 
to  the  free  vertex  of  T+.  Similar  remarks  apply  to  points  in  T~  and  the  position  vector  p~  , 
except  that  the  latter  is  directed  toward  the  free  vertex  of  T~ .  The  plus  or  minus  designation  of 
a  tetrahedron  is  determined  by  the  choice  of  a  positive  flux  reference  direction,  which  is  assumed 
to  be  from  T+  to  T~ .  The  basis  functions  that  are  used  for  volume  part  of  the  moment  method 
solution  of  (2.15)  are  such  that  one  basis  function  is  associated  with  each  face  of  the  tetrahedral 
model  of  V.  The  basis  function  associated  with  the  nth  face  is 


fn(r ) 


a"  n+ 
3V+  Pn  ’ 

r  E  T+- 

'  c  ±n  i 

-Thx-  p~ 
3Vn  1  n  ’ 

r  E 

0, 

otherwise. 

(A.21) 


where  an  is  the  area  of  the  face  and  V±  is  the  volume  of  T±.  (Note  the  convention  that  subscripts 
refer  to  faces  while  superscripts  refer  to  tetrahedrons.) 


The  unknown  volume  electric  flux  density  that  enters  (2.15)  through  (2.11,2.13,2.14)  can  be 
represented  throughout  V  by 


N 

D(r)  =  Y/Dnfn(r) 

n=  1 


(A.22) 


The  summation  is  over  the  N  faces  that  make  up  the  tetrahedral  model  of  V. 


The  basis  functions  fn(r )  have  several  properties  that  make  them  useful  for  representing  D[r). 

1)  Within  each  tetrahedron  D(r)  is  the  sum  of  four  linearly  independent  basis  functions  (one 
associated  with  each  face)  which  can  be  combined  to  represent  a  constant  vector  in  any  direction. 
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2)  fn(r)  has  no  component  normal  to  any  face  except  the  common  face  of  the  conjoined  pair 

Tn  and  Tn  ■ 

3)  The  component  of  fn  normal  to  the  nth  face  is  constant  and  continuous  across  the  face 
because  the  normal  component  of  p±  along  face  n  is  just  the  height  of  with  face  n  as  the  base 
and  the  height  expressed  as  3V±/an.  This  latter  factor  normalizes  fn  in  (A. 21)  such  that  its  flux 
density  normal  to  face  n  is  unity,  ensuring  continuity  of  the  component  of  fn,  normal  to  the  face. 


4)  The  divergence  of  the  basis  function  is 


V  •  fn(r)  =  { 

where  the  divergence  in  is  given  by 


an  r  rz  rT'Jr 

\/+  5  '  ^  -1  n  ’ 

v  n 

_  an  r  rz  rT'~ 

j  '  ^  n  > 

vn 

0,  otherwise. 


(P 


r2  diHPn)2fn 
dpn 


(A.23) 


(A.24) 


V  •  J(r)  =  —jcop(r)  (A.25) 

The  charge  density  is  proportional  to  V  •  fn  through  (A.25),  (2.16)  and  (A.22)is  constant  within 
each  tetrahedron. 

5)  The  moment  of  f  n  over  is 


/  fndv  =  (A-26) 

where  pcp^  is  the  vector  between  the  free  vertex  and  the  centroid  of  with  pc~  directed  toward 
and  pjj+  directed  away  from  the  vertex.  (A. 26)  may  be  most  easily  derived  by  expressing  the 
integral  in  terms  of  volume  coordinates. 

6)  If  face  n  is  on  the  boundary  of  V,  then  only  one  of  the  tetrahedrons,  T+  or  T~,  is  interior 
to  V.  In  this  case  it  is  assumed  that  fn,  is  defined  only  over  the  interior  tetrahedron  and  that  the 
exterior  tetrahedron  is  not  defined. 

Some  of  the  less  obvious  of  these  properties  can  be  verified.  An  important  interpretation  of 
(A. 22)  that  follows  from  1)  and  3)  is  that  the  expansion  coefficient  Dn,  represents  the  normal 
component  of  D{r )  at  the  nth  face. 

The  modeling  properties  of  the  basis  functions  f  n  are  similar  to  those  of  the  rooftop  functions 
that  are  used  to  model  surface  current  distributions.  Property  1)  guarantees  that  the  expansion 
(A. 22)  can  provide  at  least  a  piecewise  constant  approximation  to  the  flux  density.  Furthermore, 
each  function  varies  linearly  with  distance  from  its  defining  vertex.  This  provides  some  capability  to 
model  linear  field  variations.  However,  like  the  rooftop  functions,  these  basis  functions  are  constant 
in  the  transverse  direction  so  they  cannot  accurately  represent  a  vector  field  with  arbitrary  linear 
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variation;  they  can  reproduce  linear  variations  of  the  field  only  if  the  variation  is  parallel  to  the 
direction  of  the  field  defined  by  the  basis  function.  In  order  to  solve  (2.15)  with  D(r)  represented 
by  (A. 2),  it  is  necessary  to  find  expressions  for  «7(r)  and  p(r)  in  terms  of  the  basis  functions  and 
expansion  coefficients.  From  (A. 22)  and  (2.16) 

N 

J{r)  =  ju]^DnK{r)fn(r)  (A.27) 

n= 1 

In(A.27),  the  parameter  n(r)  is  taken  inside  the  summation  to  emphasize  that  the  polarization  or 
contrast  current  at  each  point  in  the  body  is  a  weighted  sum  of  the  basis  functions  at  that  point 
times  a  factor  that  describes  the  medium  at  that  point.  By  substituting  (A.27)  into  (A. 25)  and 
using  (A. 23),  the  charge  density  is  found  to  be  represented  by 


p{r) 


N  N 

Y  DnK(r)V  •  fn(r)  -  Y  Dnfn(r )  •  V/e(r) 


n= 1 


n= 1 


(A.28) 


The  first  summation  is  the  induced  volume  charge  densities  pvn  associated  with  the  basis  functions 
and  can  be  evaluated  with  the  aid  of  (A. 23) 


Pvn(r) 


_  r)  K+  an  r  fz  rT'~ 

LJn'^jn  ,r+  j  '  ^  ±n 

vn 

<  Dnn~f^,  r  £  T~ 

v  n 

0,  otherwise 


(A.29) 


where  is  the  constant  value  of  k(t)  in  T±.  The  second  summation  in  (A.28)  represents  the 
induced  surface  charge  densities  associated  with  the  basis  functions.  When  the  nth  face  separates 
dissimilar  media,  k  is  discontinuous  and  its  gradient  is  a  generalized  function  representing  a  surface 
charge  density  psn,  residing  on  an, 


Psn(r ) 


Du(k+  -  Kn),  r  E  an 
0,  otherwise 


(A.30) 


An  induced  volume  charge  associated  with  each  basis  function  exists  in  each  tetrahedron  for  which 
the  medium  parameters  differ  from  those  of  free  space.  It  can  be  shown  that  the  formulation  of  the 
problem  in  (2.11,  2.13,  2.14,  2.15,  A. 25)  guarantees  that  the  total  volume  charge  in  a  homogeneous 
region  is  zero.  In  the  numerical  procedure,  the  total  volume  charge,  which  is  the  sum  of  the  volume 
charges  associated  with  the  four  basis  functions  in  each  tetrahedron,  is  not  explicitly  constrained  to 
be  zero,  but  should  be  small  for  any  solution  that  is  close  to  the  true  solution.  An  induced  surface 
charge  exists  only  on  faces  that  separate  dissimilar  media,  and  the  charge  density  is  constant  on 
such  faces.  The  correct  surface  charge  representation  is  obtained  at  the  faces  because  of  the  careful 
construction  of  the  representation  of  J.  Equations  (A. 22),  (A.27),  and  (A.28)  can  be  used  in 
(2.11)-(2.14)  to  obtain  an  equation  involving  only  the  unknown  coefficients  Dn  and  known  vector 
functions. 


A. 2. 2  Testing  Procedure  of  SWG  Basis  Function 

The  next  step  in  applying  the  method  of  moments  is  to  select  a  testing  procedure  that  will  generate 
N  independent  equations  for  the  unknown  expansion  coefficients.  The  expansion  functions  fm 
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developed  in  the  previous  section  are  chosen  as  testing  functions.  With  a  symmetric  product 
defined  as 

<  f,  9  >=  f  f-gdv  (A.31) 

Jv 

(2.15)  is  tested  with  fm,  yielding 

<  D/e,  f  rn  >  +jui  <  A,  fm  >  +  <  V4>,  fm  >=<  Emc,  fm  >,  m  =  1, 2,  •  •  •  ,  N.  (A.32) 

This  represents  N  equations  for  the  N  unknown  coefficients  {Dn}  and  can  be  written  in  matrix 
form  as 


[Smn\[Dn\  =  [Em\  (A.33) 

where  [Smn]  is  an  N  x  N  matrix  with  dimensions  of  m4/F,  and  [Dn\  and  [Em]  are  column  vectors 
of  length  N.  Since  the  elements  of  [Smn]  and  [Emn]  cannot  be  evaluated  in  closed  form,  it  is 
necessary  to  derive  accurate  and  efficient  numerical  approximations  for  these  quantities. Formulas 
for  the  matrix  (A.33)  elements  and  some  considerations  relevant  to  their  computation  are  given  in 
the  next  section 


A. 2. 3  Elements  of  the  Matrix  and  Excitation  Vector  of  SWG 


The  formulas  that  are  used  to  fill  the  matrix  [5mn]  and  the  excitation  vector  [Em\  are  derived  from 
(A.32)  as  follows.  The  first  term  of  (A.32)  can  be  written  as 
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< 


D/i,  fm  >  =  £  A*  <  fje  > 


71=1 

N 


£a*  4  / 


1 


fn  ■  f  mdv  +  —  /  f  n  '  fmdv 

£ll  JTry, 


(A.34) 


(Recall  that  fm,  is  zero  outside  of  These  integrals  can  be  evaluated  by  utilizing  normalized 
volume  coordinates.  By  using  the  volume  coordinate  representation  of  the  basis  functions  and 
formulas  for  integration  in  volume  coordinates,  the  integral  of  fifj  over  Tv  can  be  evaluated  as 


JTp  fi  ■  fjdv  =  ^[^|**cp|2  -  (r<  +  rj)  ■  rcp  +  (n  ■  rj) 

-^(ri  ■  r2  +  ri  ■  r3  +  ri  ■  r4  +  r2  ■  r3  +  r2  ■  r4  +  r3  ■  r4)]  (A.35) 

where  rcpis  the  vector  from  O  to  the  centroid  of  Tp  and  r i  ,  r2,  r3,  r4  are  vectors  from  O  to  the 
vertices  of  T.  This  result  can  be  used  to  complete  the  evaluation  of  (A.34).  The  second  term  and 
the  right  side  of  (A.32)  both  involve  an  integration  over  a  vector  field,  and  they  can  be  treated 
similarly.  In  particular,  A  and  Etnc  are  assumed  to  vary  slowly  enough  to  be  approximated  within 
each  tetrahedron  by  their  values  at  the  centroid  of  the  tetrahedron.  Then,  by  (A. 26), 


/  (  A  \  ,  \  ^  f  A(r^f)  \  dm  Pm  I  (  J4.(r‘m.  )  1  amPm 
\V  Einc  J'Jm/  ~  V  Einc{rc+)  )  '  3  V  EinC{rCm)  )  '  3 


(A.36) 
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If  T+  or  Tm  is  not  in  V,  then  the  corresponding  term  is  absent  in  (A. 36).  In  (A. 36),  is  the 

vector  from  O  to  the  centroid  of  and  p^  is  the  vector  from  the  vertex  opposite  am  in  Tm  to 
the  centroid  of  Tm. 


The  scalar  potential  term  of  (A.  32)  can  be  written  as 

<v<^/m>=  (  f  <S>V-fmdv  (A.37) 

Js  Jv 

where  Sis  the  boundary  of  V.  Since  the  tangential  component  of  fm,  is  discontinuous  at  tetrahedral 
boundaries,  verification  of  (A.37)  requires  that  the  contributions  from  individual  tetrahedral  volume 
elements  and  their  boundaries  be  summed.  This  leads  directly  to  the  last  term  of  (A.37).  The 
surface  integral  follows  from  the  fact  that  the  normal  component  of  f  m ,  is  continuous  or  vanishes 
at  all  boundaries  internal  to  V.  Then,  with  the  assumption  that  $  is  sufficiently  slowly  varying 
that  it  may  be  replaced  by  its  value  at  the  centroid  of  a  tetrahedron  or  face, (A.37)  becomes 


f  am[^(rcm  )  -  $«+)],  T+  and  Tm  in  V 

1  ±am[$(r*m  ±  1/3 p^)  -  $(r£±)],  T^not  inC 


(A.38) 


The  argument  (r^f  ±  1/3 pc±)  denotes  that  is  evaluated  at  the  centroid  of  face  m  when  is 
not  in  V.  When  the  elements  of  [Smn]  are  calculated,  the  contributions  to  A  and  $  from  a  single 
basis  function  are  needed.  These  are  given  by 
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(A.40) 


where  is  taken  to  be  zero  if  is  not  in  V. 


By  expressing  them  in  terms  of  normalized  volume  coordinates,  the  integrals  in  (A.39)  and 
(A.40)  may  be  readily  evaluated  numerically.  They  can,  in  fact,  be  expressed  in  terms  of  four 
independent  scalar  integrals  over  each  tetrahedron.  These  four  integrals  contribute  to  16  different 
elements  of  [*S'mn]  •  Considerable  savings  in  computation  time  is  achieved  by  making  use  of  these 
multiple  contributions. 
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